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ABSTRACT 

We present an open-source Thermochemical Equilibrium Abundances (TEA) code that calculates the abun¬ 
dances of gaseous molecular species. The code is based on the methodology of White et al. (1958) and Eriksson 
(1971). It applies Gibbs free-energy minimization using an iterative, Lagrangian optimization scheme. Given 
elemental abundances, TEA calculates molecular abundances for a particular temperature and pressure or a 
list of temperature-pressure pairs. We tested the code against the method of Burrows & Sharp (1999), the 
free thermochemical equilibrium code CEA (Chemical Equilibrium with Applications), and the example given 
by White et al. (1958). Using their thermodynamic data, TEA reproduces their final abundances, but with 
higher precision. We also applied the TEA abundance calculations to models of several hot-Jupiter exoplan¬ 
ets, producing expected results. TEA is written in Python in a modular format. There is a start guide, a user 
manual, and a code document in addition to this theory paper. TEA is available under a reproducible-research, 
open-source license via https : / /github . com/dzesmin/TEA. 

Subject headings: astrochemistry - molecular processes - methods: numerical - planets and satellites: atmo¬ 
spheres - planets and satellites: composition - planets and satellites: gaseous planets 


1. INTRODUCTION 

There are two methods to calculate chemical equilibrium: 
using equilibrium constants and reaction rates, i.e., kinetics, 
or minimizing the free energy of a system (Bahn & Zukoski 
1960, Zeleznik & Gordon 1968). 

The kinetic approach, where the pathway to equilibrium 
needs to be determined, is applicable for a wide range of tem¬ 
peratures and pressures. However, using kinetics for thermo¬ 
chemical equilibrium calculations can be challenging. Chem¬ 
ical equilibrium can be calculated almost trivially for several 
reactions present in the system, but as the number of reactions 
increases, the set of numerous equilibrium constant relations 
becomes hard to solve simultaneously. To have an accurate ki¬ 
netic assessment of the system, one must collect a large num¬ 
ber of reactions and associate them with the corresponding 
rates. This is not an issue at lower temperatures, where re¬ 
action rates are well known. However, at high temperatures, 
where thermochemical equilibrium should prevail, one needs 
to know forward and reverse reactions and corresponding re¬ 
action rates, which are less well known. 

The advantage of the free energy minimization method is 
that each species present in the system can be treated inde¬ 
pendently without specifying complicated sets of reactions a 
priori, and therefore, a limited set of equations needs to be 
solved (Zeleznik & Gordon 1960). In addition, the method 
requires only knowledge of the free energies of the system, 
which are well known, tabulated, and can be easily interpo¬ 
lated or extrapolated. 

Thermochemical equilibrium calculations have been 
widely used in chemical engineering to model combus¬ 
tion, shocks, detonations and the behaviour of rockets and 
compressors (e.g.. Miller et al. 1990, Belford & Strehlow 
1969). In astrophysics, they have been used to model the 
solar nebula, the atmospheres and circumstellar envelopes 
of cool stars, and the volcanic gases on Jupiter’s satellite lo 
(e.g., Lauretta et al. 1997, Lodders & Eegley 1993, Zolotov 
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& Eegley 1998). Thermochemistry also governs atmospheric 
composition in vast variety of giant planets, brown dwarfs, 
and low-mass dwarf stars (Lodders & Eegley 2002, Visscher 
et al. 2006, 2010a, and references therein ). 

1.1. Chemical Models of Exoplanets 

To perform a comprehensive study of a planetary atmo¬ 
sphere, aside from thermoequilibrium chemistry, one must 
consider disequilibrium processes like photochemistry, ver¬ 
tical mixing, horizontal transport, and transport-induced 
quenching (Moses et al. 2011, 2013, Venot et al. 2012, 2014, 
Agtindez et al. 2012, 2014, Line et al. 2010, 2011, Visscher 
& Moses 2011). Today, we have ID chemical models that in¬ 
tegrate thermochemistry, kinetics, vertical mixing, and photo¬ 
chemistry (Line et al. 2011, Moses et al. 2011, Visscher et al. 
2010b). These models have an ability to smoothly transition 
from the thermochemical-equilibrium regime to transport- 
quenched and photochemical regimes. Specifically, in giant 
planets, we can distinguish three chemical layers: deep within 
the planetary atmosphere, the temperatures and pressures are 
so high that chemical reaction timescales are short, ensuring 
a chemical equilibrium composition; at lower temperatures 
and pressures higher in the atmosphere, the timescales for 
chemical reactions slows down, reaching the vertical trans¬ 
port timescale and smoothing the vertical mixing-ratio profile 
by producing quenched abundances; high in the atmosphere, 
the host star’s ultraviolet radiation destroys stable molecules, 
driving photochemical reactions. 

Photochemical models today face several difficulties. They 
lack high-temperature photochemical data, and the list of re¬ 
actions and associated rate coefficients are not well defined 
or are conflicted (Venot et al. 2012, Visscher et al. 2010b). 
In addition, the exoplanet photospheres observed with current 
instruments are sampled within the region of the atmosphere 
dominated by vertical mixing and quenching, but not by pho¬ 
tochemistry (Line & Yung 2013). 

The majority of early hot-Jupiter atmospheric models as¬ 
sumed chemical composition consistent with thermochemical 
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equilibrium (e.g., Burrows et al. 2007, Fortney et al. 2005, 
Marley et al. 2007, Fortney et al. 2010, Burrows & Sharp 
1999, Sharp & Burrows 2007, Rogers et al. 2009). More re¬ 
cently, a variety of ID photochemical models has been used to 
explore the compositions of hot Jupiters (Moses et al. 2011, 
Zahnle et al. 2009, Line et al. 2010, Kopparapu et al. 2012, 
Venot et al. 2012, Line & Yung 2013, Visscher et al. 2006, 
2010b, Venot et al. 2012). A common conclusion of these 
studies is that in hot atmospheres {T > 1200 K), disequilib¬ 
rium effects are so reduced that thermochemical equilibrium 
prevails. 

Using secondary eclipse observations as the most fruit¬ 
ful technique today to assess atmospheric composition (e.g., 
Knutson et al. 2009a,b, Machalek et al. 2008, Stevenson et al. 
2012, Fraine et al. 2013, Crossfield et al. 2012, Todorov 
et al. 2012, Desert et al. 2009, Deming et al. 2011, Demory 
et al. 2007, Madhusudhan & Seager 2009, Madhusudhan et al. 
2011a, Blecic et al. 2013, 2014), Line & Yung (2013) studied 
the most spectroscopically active species in the infrared on 
eight hot planets (GJ436b, WASP-12b, WASP-19b, WASP- 
43b, TrES-2b, TrES-3b, HD 189733b, and HD 149026b), 
with equilibrium temperatures ranging between 744 K and 
2418 K. They chose to evaluate the presence of disequilibrium 
chemistry at 100 mbar, where most secondary-eclipse obser¬ 
vations sample, i.e., where their thermal emission weighting 
functions usually peak. They find that all of the models are 
consistent with thermochemical equilibrium within 3 (t (the 
work of Stevenson et al. 2010, however, questions this conclu¬ 
sion for GJ436b). They also show that for the hottest planets, 
(Tioomb > 1200 K), CH 4 , CO, HjO, and Hj should be in ther¬ 
mochemical equilibrium even under a wide range of vertical 
mixing strengths. 

Thermochemical equilibrium calculations are the starting 
point for initializing models of any planetary atmosphere. 
In general, thermochemical equilibrium governs the compo¬ 
sition of the deep atmospheres of giant planets and brown 
dwarfs, however, in cooler atmospheres thermoequilibrium 
calculations are the necessary baseline for further disequilib¬ 
rium assessment. They can also provide a first-order approxi¬ 
mation for species abundances as a function of pressure, tem¬ 
perature, and metallicity for a variety of atmospheres (e.g., 
Visscher et al. 2010b, Lodders & Eegley 2002). 

The Gibbs free energy minimization method for calculat¬ 
ing thermochemical equilibrium abundances of complex mix¬ 
tures was first introduced by White et al. (1958). Prior to 1958 
all equilibrium calculations were done using equilibrium con¬ 
stants of the governing reactions. White et al. (1958) were 
the first to develop a method that makes no distinction among 
the constituent species and does not need a list of all possible 
chemical reactions and their rates. Rather, it depends only on 
the chemical potentials of the species involved. To derive the 
numerical solution, they apply two computational techniques: 
a steepest-descent method applied to a quadratic fit and the 
linear programming method. 

Eollowing their methodology, Eriksson (1971) developed 
the SOLGAS code that calculates equilibrium composition 
in systems containing ideal gaseous species and pure con¬ 
densed phases. Subsequent modification of this code were 
made by Eriksson (& Rosen (1973), Eriksson (1975), and Be- 
smann (1977), after which the code was modified for astro- 
physical applications and called SOLAGASMIX by Sharp & 
Huebner (1990), Petaev & Wood (1998), Burrows & Sharp 
(1999), and Sharp & Burrows (2007). 

The Gibbs free energy minimization approach has been 


used by many authors in the exoplanetary field (e.g., Sea¬ 
ger 1999, 2010, Madhusudhan & Seager 2009, 2010, Sharp 
& Huebner 1990). In addition to SOLAGASMIX and 
other proprietary codes (e.g., CONDOR by Eegley & Lod¬ 
ders 1994, Lodders & Eegley 1994), and one analytic 
method to calculate major gaseous species in planetary at¬ 
mospheres by Burrows & Sharp (1999), just one free- 
software code CEA, (Chemical Equilibrium with Applica¬ 
tions, http://WWW.grc.nasa.gov/WWW/CEAWeb, by 
Gordon & McBride 1994), is available to the exoplanet com¬ 
munity. 

In this paper, we present a new open-source 
code. Thermochemical Equilibrium Abundances 
(TEA). The TEA code is a part of the open-source 
Bayesian Atmospheric Radiative Transfer project 
(https : //github . com/ joeharr4 /BART). This 
project consists of three major parts: TEA - this code, a 
radiative-transfer code that models planetary spectra, and 
a statistical module that compares theoretical models with 
observations. 

TEA calculates the equilibrium abundances of gaseous 
molecular species. Given a single T,P point or a list of 
T, P pairs (the thermal profile of an atmosphere) and elemen¬ 
tal abundances, TEA calculates mole fractions of the desired 
molecular species. The code is based on the Gibbs free energy 
minimization calculation of White et al. (1958) and Eriksson 
(1971). TEA uses 84 elemental species and the thermody¬ 
namical data for more then 600 gaseous molecular species 
available in the provided JANAE (Joint Army Navy Air Eorce) 
tables (http : //kinetics . nist. gov/ janaf /, Chase 
et al. 1982, Chase 1986). TEA can adopt any initial elemental 
abundances. Eor user convenience a table with solar photo- 
spheric elemental abundances from Asplund et al. (2009) is 
provided. 

The code is written in Python in an architecturally modular 
format. It is accompanied by detailed documentation, a 
start guide, the TEA User Manual (Bowman and Ble¬ 
cic), the TEA Code Description document (Blecic and 
Bowman), and the TEA Theory document (this paper), 
so the user can easily modify it. The code is actively 
maintained and available to the scientific community via 
the open-source development website GitHub. com 
(https://github.com/dzesmin/TEA, 
https://github.com/dzesmin/TEA-Examples). 
This paper covers an initial work on thermochemical cal¬ 
culations of species in gaseous phases. Implementation of 
condensates is left for future work. 

In this paper, we discuss the theoretical basis for the method 
applied in the code. Section 2 explains the Gibbs Eree en¬ 
ergy minimization method; Section 3 describes the general 
Lagrangian optimization method and its application in TEA; 
in Section 4 we introduce the Lambda Correction algorithm 
for handling negative abundances that follow from the La¬ 
grangian method; Section 5 describes the layout of the TEA 
code; Section 6 explores chemical equilibrium abundance 
profiles of several exoplanetary atmospheres; Section 7 com¬ 
pares our code to other methods available, and Section 9 states 
our conclusions. 

2 . GIBBS EREE ENERGY MINIMIZATION METHOD 

Equilibrium abundances can be obtained by using different 
combinations of thermodynamical state functions: tempera¬ 
ture and pressure - (f,p), enthalpy and pressure - en¬ 

tropy and pressure - {S,p), temperature and volume - (f,v). 
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internal energy and volume - (t/,v), etc. Depending on how 
the system is described, the condition for equilibrium can be 
stated in terms of Gibbs free energy, helmholtz energy, or en¬ 
tropy. If a thermodynamic state is defined with temperature 
and pressure, Gibbs free energy (G) is most easily minimized, 
since those two states are its natural, dependent variables. 

Gibbs free energy represents a thermodynamic potential 
that measures the useful work obtainable by the system at a 
constant temperature and pressure. Thus, the Gibbs free en¬ 
ergy minimization method minimizes the total chemical po¬ 
tential of all involved species when the system reaches equi¬ 
librium. 

The Gibbs free energy of the system at a certain temperature 
is the sum of the Gibbs free energies of its constituents: 


GUT) 

RT 




r g°(7^) 

. RT 


-l-lnP-l-ln- . 

Ni 


(9) 


Equation (9) requires a knowledge of the free 
energy of each species as a function of tempera¬ 
ture. These can be obtained from the JANAE tables 
(http : //kinetics . nist. gov/ janaf/, Chase et al. 
1982, Chase 1986, Burrows & Sharp 1999), or easily derived 
from other tabulated functions. 

To extract free energies, g^{T)/RT, from the JANAE tables, 
we used the expression given in Eriksson (1971), Equation 
( 2 ): 


GUT) = Y.^‘(T), (1) 

i 

where Gsys{T) is the total Gibbs free energy of the system 
for n chemical species, Gi{T) is the Gibbs free energy of a 
gas species i, and T is the temperature. The total Gibbs free 
energy of the system is expressed as the sum of the number 
of moles X of the species /, x,, and their chemical potentials 
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where g°(7’) is given in J/mol, R = 8.3144621 J/K/mol,is 
the enthalpy (heat content) in the thermodynamical standard 
state at a reference temperature of 25°C = 298.15 K, G? is the 
Gibbs free energy in J/mol, (G^-Z/jgg/T) is the free-energy 
function in J/K/mol, and AyZ/jgg is the heat of formation at 
298.15 K in kJ/mol. Thus, our conversion equation becomes: 


n 

GUT) = Ylxigi(T). ( 2 ) 

i 

The chemical potential gi{T) depends on the chemical poten¬ 
tial at the standard state g^{T) and the activity a,, 

gi(T)=g^(T) + RT\nai, (3) 

where R is the gas constant, R = ksNA, and kg and Na are 
the Boltzmann constant and Avogadro’s number, respectively. 
Activities for gaseous species, which are treated as ideal, are 
equal to the partial pressures, and for condensates they equal 
1 : 


fl, =Pi = P^, forgases 

(4) 

a, = 1 , forcondensates, 

(5) 

where P is the total pressure of the atmosphere, N is the total 
number of moles of all species involved in the system. Hence, 
Equation (3) for gaseous species becomes: 

gi(T) = g°(T)+RT\nPi. 

( 6 ) 

Combining Equation ( 6 ) with Equation (2), the Gibbs 
ergy of the system becomes: 

free en- 

n 

GUT) = Y,Xi(g^U) + RT\nP^ , 

(7) 


or. 


G.,.(7’) = ^x,(g°(r) + J?7’lnP+/?7’ln^) , (8) 

i 

Eor our purposes, it is more convenient to write Equation ( 8 ) 
in unitless terms: 
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G|’-// 29 g/T is the fourth term in the JANAE tables and 
A/HU is the sixth. The free energy function of a species 
corresponding to a temperature other than those provided in 
the JANAE tables is calculated using spline interpolation. 

Alternatively, the free energies can be calculated using the 
eighth term in the JANAE tables, following Equation (3) from 
Eriksson (1971): 


In (10) log UKf), ( 12 ) 

RT 

where Kf is the equilibrium constant of formation. 

To determine the equilibrium composition, we need to find 
a non-negative set of values x, that minimizes Equation (9) 
and satisfies the mass balance constraint: 


'^aijXi = bj, (j=l,2,...,m), (13) 

;=1 

where the stoichiometric coefficient indicates the number 
of atoms of element j in species i (e.g., for CH 4 the stoichio¬ 
metric coefficient of C is 1 and the stoichiometric coefficient 
of H is 4), and bj is the total number of moles of element j 
originally present in the mixture. 

We use the reference table containing elemental solar abun¬ 
dances given in Asplund et al. (2009) Table 1 for b values. 
Asplund et al. (2009) adopt the customary astronomical scale 
for logarithmic abundances, where hydrogen is defined as log 
eH = 12.00, and log eX = \og{Nx /Nh)+ 12, where Nx and Nh 
are the number densities of element X and H, respectively. 
Thus, their values are given in dex (decimal exponent) units. 
We transform these values into elemental fractions by num¬ 
ber, i.e., ratio of number densities. We convert each species 
dex elemental abundance into number density and divide it by 
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the hydrogen number density (Asplund et al. 2009, Section 3). 
The final output are fractional abundances (mole mixing frac¬ 
tions), i.e., the ratio of each species’ number of moles to the 
number of moles in the mixture. 

3. LAGRANGIAN METHOD OF STEEPEST DESCENT 

To find equilibrium abundances of the desired molecu¬ 
lar species at a given temperature and pressure, we need 
to minimize Equation (9). To do so, we have to apply a 
technique that minimizes a multi-variate function under con¬ 
straint. There are many optimization techniques used to find 
the minima of a function subject to equality constraints (e.g., 
line search method, Dantzig-simplex method for linear pro¬ 
gramming, Newton-Raphson method, Hessian-conjugate gra¬ 
dient method, Lagrangian steepest-descent method). The 
main advantage of the Lagrangian steepest-descent method is 
that the number of equations to solve scales with the number 
of different types of atoms present in the mixture, which is 
usually a much smaller number than the possible number of 
molecular constituents. This allows the code to be executed 
much faster than in other methods. 

Gradient descent, also known as steepest descent, is an al¬ 
gorithm for finding a local minimum of a function. At each 
iteration, the method takes steps towards the minimum, where 
each step is proportional to the negative gradient of the func¬ 
tion at the current point. If a function f{x) is defined and 
differentiable in the neighborhood of a point a, then /(x) de¬ 
creases most rapidly in the direction of the negative gradient, 
-V/(a). From this, it follows that if b = a- AV/(a), then 
f{a) > f{b) if A is small enough. Starting with a guess xq for 
a local minimum of /, and considering a sequence xo, xi, X 2 ,... 
such that x„+i = x„ - AV/(x„), n > 0 , one gets /(xq) > /(xi) 
> fix 2 ) >... . This sequence of x„ converges to a desired lo¬ 
cal minimum if the correct A value is assigned. The value of 
A can vary at each iteration. If the function / is convex, the 
local minimum is also the global minimum. 

Our code implements a more complex version of the 
method outlined above. The problem consists of some func¬ 
tion /(x,y) subject to a constraint g(x,y) = C. In this case, we 
need both / and g to have continuous first partial derivatives. 
Thus, we introduce a new variable called the Lagrangian mul¬ 
tiplier, TT, where; 


A(x, y, A) = /(x, y) ± TT (g (x, y) - C), (14) 

which allows us to find where the contour of g{x,y) = C tan¬ 
gentially touches fix,y) (Figure 1). The point of contact is 
where their gradients are parallel: 

yxyf{x,y) = -nVxyg(x,y). (15) 

The constant tt allows these gradients to have different mag¬ 
nitudes. To find the minimum, we need to calculate all partial 
derivatives of the function A, equate them with zero, 

Vj:,j,,^A(x,y,7r) = 0 , (16) 

and follow the same iteration procedure as explained above. 

3.1. Lagrangian Method in TEA 

To implement this in our code, we followed the methodol¬ 
ogy derived in White et al. (1958). We applied an iterative 
solution to the energy minimization problem, where the mole 


numbers of the desired molecular species are recomputed at 
each step and the new direction of steepest descent is calcu¬ 
lated. This produces improved mole number values, which 
however, could be negative. Thus, two short procedures are 
required in each iteration cycle: solving a set of simultane¬ 
ous linear equations for an improved direction of descent (de¬ 
scribed in this Section) and approximately minimizing a con¬ 
vex function of one variable. A, to ensure that all improved 
mole number values are positive (Section 4). 

To calculate the direction of steepest descent (following the 
methodology derived in Section 3) and initiate the first itera¬ 
tion cycle, we first need to solve the mass balance equation. 
Equation (13). We start from any positive set of values for the 
initial mole numbers, y = (yi ,y 2 , ...,y„), as our initial guess; 


'^aijyi = bj (;'= l,2,...,m). (17) 

1=1 

To satisfy the mass balance Equation (17), some y, variables 
must remain as free parameters. In solving these equations, 
we leave as many free parameters as we have elements in the 
system, thus ensuring that the mass balance equation can be 
solved for any number of input elements and output species 
the user chooses. We set all other y, to a known, arbitrary 
number. Initially, the starting values for the known species 
are set to 0.1 moles, and the mass balance equation is calcu¬ 
lated. If that does not produce all positive mole numbers, the 
code automatically sets known parameters to 10 times smaller 
and tries again. The initial iteration input is set when all mole 
numbers are positive, and the mass balance equation is satis¬ 
fied. 

To follow with the Lagrangian method, we denote two 
terms in Equation (9) as; 


/(T) 

c,= ^^-l-lnP, (18) 

RT ’ 

where P is the pressure in bar. Using a, we denote the right 
side of Equation (9) as the variable fi{Y): 


fi(Y) = yi 



(19) 


where Y = (yi,y 2 ,...,y„) and y is the total initial number of 
moles. The left side of Equation (9), Gsys(T)/RT, we denote 
as function F{Y): 


TO = ^y/[Q + ln^ 
1=1 ^ 


( 20 ) 


Then, we do a Taylor series expansion of the function F 
about Y. This yields a quadratic approximation QiX): 


Q{X) = F{X) 


X=Y 


■E 


dF 

dxi 


X=Y 




lEi: 


d^F 

dxidxk 


AiAk. 

X=Y 


( 21 ) 


where A, = Xi-yi, and x; are the improved mole numbers. 
This function is minimized using the Lagrangian principle. 
We now introduce Lagrangian multipliers as tt^: 
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Fig. 1.— Example of the Lagrangian minimization approach. Left: The 3D illustration of the minimization problem. Blue lines indicate starting and ending 
values of f{x,y) during minimization. An (x, y) pair is found that minimizes f(x,y) (bottom blue line) subject to a constraint g(x,y) = C (red line). Right: Contour 
map of the left figure. The point where the red line (constraint) tangentially touches a blue contour is the solution. Since di < do, the solution is the minimum of 


GiX) = Q{X) + Y, T^ji- Y ’ ( 22 ) 

j ‘ 

and calculate the first derivatives, dGjdxi, of the new func¬ 
tion. We equate them to zero to find the minima, dGjdxi = 
0 . 

We solve for x,- from Equation (22) by combining Equation 
(17) and (19) with the fact that x is the sum of the improved 
mole numbers, x = The improved number of moles, 

X/, are given as: 


Xi = -fi(Y) + {^)x+{Y'^j^ij)yi’ ( 23 ) 

^ 1=1 

while the Lagrangian multipliers, nj, are expressed as: 


n 

rnT^\_ + ri2'K2 + - + r\mT^m -I-m = ^ fl/i /;(T), 

/=! 

n 

r2\ TTi -I- r227r2 -I-... -I- r2m7r„, -I- ^2 « = ^ aa fi(Y ), 

/=! 

■ 5 

, (27) 

• 5 

n 

^nil'll mm'^m ^ ^ ^ ^im fiiX^ ^ 

i=l 

n 

/?! TTi +/727r2 + • ■ • + + 0 M ^/;(y) , 

(=1 

where: 


m n 




n 

^ijYi = E^' 

i=l 


r g°(r) 
. RT 


-l-lnP-l-ln^j , 


(24) 


where j iterates over the m elements and i iterates over the n 
species, x and y are the sums of improved and initial num¬ 
ber of moles, respectively. Using Equation (19), we can now 
rewrite Equation (24) as: 


T.^jbj = Yf‘(Y'^- (25) 

1=1 1=1 

If we further denote the constants with: 


n 

r jk = rkj = Y^'^'j '^'*) Y' ’ (2^) 

(=1 

combining Equations (23), (25), and (26), we get the follow¬ 
ing system of m-l-1 equations that can easily be solved: 


M = —1-l-x/y. (28) 

The solutions to Equations (27) and (28) will give ttj and u, 
and from them using Equation (23) we can calculate the next 
set of improved mole numbers, i.e., an improved direction of 
descent. A,- = x,-y,'. 

4. LAMBDA CORRECTION ALGORITHM 

Solving a system of linear equations (i.e., performing the 
Lagrangian calculation) can also lead to negative mole num¬ 
bers for some species, so a short additional step is needed to 
eliminate this possibility and guarantee a valid result. 

To do so, the difference between the initial and final values 
given by the Lagrangian calculation. A, = x, -y,, we will call 
the total distance for each species. To ensure that all improved 
mole numbers are positive, we introduce a new value. A, that 
defines the fraction of the total distance as AA, (see Eigure 2). 

The computed changes, AA,, are considered to be direc¬ 
tional numbers indicating the preferred direction of descent 
the system moves to. Other than providing all positive mole 
numbers, we determine the value A so that the Gibbs energy 
of the system must decrease, i.e., the minimum point is not 
passed (see Equation 34). 

At each Lagrangian iteration cycle we start with the initial 
positive values, y, and we get the next set of improved values 
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Every new iteration starts with a different set of y,, thus 
changing the convex function F{X), Equation 33, and produc¬ 
ing a new minimum. This yields to a new A value. A will 
be found to approach unity after some number of iterations. 
Unity in A indicates the solution is near. 

We repeat the Lagrangian method and the lambda correc¬ 
tion until a pre-defined maximum number of iterations is 
met. The final abundances are given as fractional abundances 
(mole mixing fractions), i.e., the ratio of each species’ mole 
numbers to the total sum of mole numbers of all species in the 
mixture. 


Fig. 2.— Simplified illustration of the lambda correction algorithm. Initial 
values for one hypothetical Lagrangian iteration cycle, y,, are given in green. 
These values are all positive and satisfy the mass balance equation. Equation 
17. The x; values, given in blue, are the values produced by the Lagrangian 
calculation. These values can be negative, but they also satisfy the mass 
balance equation. The x- values, given in red, are produced by choosing the 

maximum value of lambda that ensures all positive and non-zero x-. These 
values become the new initial values of y, for the next iteration cycle. 

Xi given as: 


Xi=yi + Ai. (29) 

Since we do not want any x; to be negative, the variable A 
performs a small correction: 

Xi=yi + XAi. (30) 

A takes values between 0 and 1, where value of zero implies 
no step is taken from the iteration’s original input, y,, and one 
implies that the full Lagrangian distance is travelled. A,-. We 
now rewrite Equation (19) using Equation (30) as: 

/i(A') = x;(^+lnE+ln|-), (31) 

which can be written in the form: 

/;(A) = (y, + AA,)(^ + lnP + ln|i^), (32) 

where A = y-x. Summing over i, we get a new function, 

^’(A): 


^’(A) = ^(y;+AA,)( 


8°(T) 

RT 


-l-lnE-bln 


y/-b AA,\ 
y-l-AA / 


(33) 


Thus, to ensure that the new corrected values x, are all 
positive, the distance travelled will be limited to fractional 
amounts defined by AA,, using the largest possible value of 
A that satisfies the conditions: 


1. The function called the directional derivative is defined 
and exists: 


dF{\) 

dX 




r g?(r) 

. RT 


-l-lnP-l-ln 


y,' + A A; ■ 

y-bAA - 


(34) 


2. The directional derivative does not become positive (the 
minimum point is not passed). 


5. CODE STRUCTURE 

The TEA code is written entirely in Python and uses 
the Python packages NumPy (http : //numpy . org/) and 
(http: //www. scipy.org/) along with SymPy, an ex¬ 
ternal linear equation solver (http : //sympy . org/). 

The code is divided into two parts: the pre-pipeline that 
makes the thermochemical data library and stoichiometric ta¬ 
bles, and the pipeline that performs abundance calculations. 
Given elemental abundances, TEA calculates molecular abun¬ 
dances for a particular temperature and pressure or a list of 
temperature-pressure pairs. Documentation is provided in the 
TEA User Manual (Bowman and Blecic) and the TEA Code 
Description (Blecic and Bowman) that accompany the code. 
Eigure 3 shows the layout of the TEA program’s flow. Its 
modules are: 

1. prepipe.py: Runs the readJANAF .py and 

makestoich. py modules and provides their 
common setup. 

2. readJANAF.py: Extracts relevant from all avail¬ 
able NIST-JANAE Thermochemical Tables and writes 
ASCII files. 

3. makestoich.py: Reads the chemical formula to ob¬ 
tain species names and their stoichiometric coefficients 
from each JANAE file, and elemental solar abundances 
from an ASCII file based on Asplund et al. (2009) Table 
1. The code produces an output file containing species, 
stoichiometric coefficients, and abundances. 

4. runsingle.py: Runs TEA for a single T,P pair. 

5. runatm.py: Runs TEA over a pre-atmosphere file con¬ 
taining a list of T,P pairs. 

6. readatm.py: Reads the pre-atmospheric file with mul¬ 
tiple TjP pairs. 

7. makeheader.py: Combines the stoichiometric infor¬ 
mation, Gibbs free energy per species at specific tem¬ 
peratures, and the user input to create a single file 
with relevant chemical informations further used by the 
pipeline. 

8 . balance.py: Uses species and stoichiometric informa¬ 
tion to establish viable, mass-balanced, initial mole 
numbers. 

9. format.py: Auxiliary program that manages in¬ 
put/output operations in each piece of the pipeline. 

10. lagrange.py: Uses data from the most recent itera¬ 
tion’s corrected mole numbers and implements the La¬ 
grangian method for minimization. Produces output 
with raw, non-corrected mole numbers for each species 
(values are temporarily allowed to be negative). 
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Pre-pipeline modules to make thermodynamic library 



Pipeline modules 


Multiple T-P 
input atm. file 

__ —►/ runatm.py 

readatm.py 

-OR- ,-, \ 

Single T-P 
input file 

-fc/runsingle-py)-^ 

makeheader.py 


-Use next line of input 



Yes 







Minimized 
abundances 
for T-P point 

Run by \ 
*\^unatm.py^ 

^Yes-»<; 

End of^^ 
s. atm. file , 

^^No 


No 


Yes 








Final Abundances 
for single T-P 


Final Abundances 
for multiple T-P 


Legend: 


User 

input 



Program 



Fig. 3.— Layout of the TEA pre-pipeline and pipeline modules. The modules have one of three roles: scientific calculation, file or data structure support, or 
execution of the calculation programs over temperature and pressure points in an iterative manner. In addition to the modules shown, TEA has three supporting 
modules: readconfig.py, makeatm.py, and plotTEA.py. All modules are described in the text. 
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Fig. 4.— The left and middle panels show the 0-rich and C-rich temperature and pressure {T—P) profile of WASP-12b from Stevenson et al. (2014a) with C/0 
= 0.5 and C/0 = 1.2 respectively. The right panel shows the T—P profile of WASP-43b from Stevenson et al. {2014b) with solar metallicity. 




Fig. 5.— Comparison between vertical thermochemical equilibrium distributions for WASP-12b 0-rich, left panel, and WASP-12b C-rich, right panel, elemen¬ 
tal abundance profile. The inset plots are the T—P profiles shown in Figure 4, left and middle panels. 


11. lambdacorr.py: Takes non-corrected mole numbers 
and implements lambda correction to obtain only valid, 
positive numbers of moles. Output is the corrected 
mole numbers for each species. 

12. iterate.py; Driver program that repeats lagrange.py 
and lambdacorr.py until a pre-defined maximum num¬ 
ber of iterations is met. 

13. readconflg.py: Reads TEA configuration file. 

14. makeatm.py: Makes pre-atmospheric file for a multi¬ 
ple 7’,Rrun. 

15. plotTEA.py: Plots TEA output, the atmospheric file 
with final mole-fraction abundances. 

6 . APPLICATION TO HOT-JUPITER ATMOSPHERES 

In this section, we illustrate several applications of the TEA 
code. We produced molecular abundances profiles for models 
of hot-Jupiter planetary atmospheres, given their temperature- 
pressure profiles. 

The temperature and pressure (T -P) profiles adopted for 
our thermochemical calculations are shown in Eigure 4. The 
left and middle panel show the T-P profiles of WASP-12b 
from Stevenson et al. (2014a) with the C/O ratio of 0.5 and 
1.2, respectively. The right panel shows the thermal profile 
of WASP-43b from Stevenson et al. (2014b) with solar metal¬ 
licity. These profiles are chosen for their relevance to atmo¬ 
spheric conditions at secondary eclipses. 


We chose elemental-abundance profiles with C/O > 1 and 
C/O < 1 and three profiles with solar, 10 times solar, and 50 
times solar elemental abundances to show the influence of the 
C/O ratio and metallicity on the chemistry and composition of 
extrasolar giant planets. 

We adopt Asplund et al. (2009) photospheric solar abun¬ 
dances as our baseline. To change the elemental abundance 
profile, set them to a certain C/O ratio, or enhance metallic¬ 
ity, we use our Python routine, makeAbun . py. This rou¬ 
tines is the part of the BART project and it is available to the 
community via Git hub. com under an open-source licence 
(https : //github . com/ joeharri /BART). Eor differ¬ 
ent metallicities, the routine multiples the elemental abun¬ 
dances of all species except for hydrogen and helium, pre¬ 
serving the ratio of major atomic species like C, N, and O. 

We chose to run the models for all plausible, spectroscopi¬ 
cally active species in the infrared relevant for hot-Jupiter at¬ 
mospheres: Hj, CO, COj, CH 4 , H 2 O, HCN, C 2 H 2 , C 2 H 4 , N 2 , 
NH 3 , HS, and H 2 S. Our input species are: H, He, C, N, O, S. 

Eigure 5 shows results for WASP-12b. Each T-P profile 
is sampled 100 times uniformly in log-pressure space. Eigure 
6 shows the TEA runs for WASP-43b with different metal¬ 
licities. This T-P profile is sampled 90 times in uniformly 
log-pressure space. 

As expected, Eigure 5 shows that H 2 O, CH 4 , CO, CO 2 , 
C 2 H 2 , C 2 H 4 , and HCN are under the strong influence of the 
atmospheric C/O ratio in hot Jupiters (e.g., Lodders & Eegley 
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Mixing Fraction 


Fig. 6.— Thermochemical equilibrium vertical distributions for different 
metallicities of WASP-43b assuming the F—P profile in Figure 4, right panel 
(profile given in inset). Three metallicity cases with (^= 1, 10, and 50 are 
shown from the top to the bottom. 

2002, Seager et al. 2005, Fortney et al. 2005, Madhusudhan 
et al. 2011a,b, Madhusudhan 2012, Madhusudhan & Seager 
2011, Moses et al. 2013). These species are plotted in solid 
lines, while species with only small influence from the C/O 
ratio are plotted as dashed lines. 

The results also show, as expected, that CO is a major at¬ 
mospheric species on hot Jupiters for all C/O ratios and metal¬ 


licities (Figures 5 and 6 ), because CO is chemically favored 
over H 2 O. Other oxygen-bearing molecules like H 2 O and CO 2 
are more abundant when C/0< 1, while CH 4 , C 2 H 2 , and C 2 H 4 
become significant species when C/0> 1. Species like N 2 and 
NH 3 that do not contain carbon or oxygen are much less af¬ 
fected by the C/O ratio. 

H 2 O is abundant in hot-Jupiter atmospheres (e.g.. Bur¬ 
rows & Sharp 1999, Fodders & Fegley 2002, Hubeny & Bur¬ 
rows 2007, Sharp & Burrows 2007) due to the large solar 
abundances of oxygen and hydrogen. Even disequilibrium 
processes like photochemistry cannot deplete its abundance. 
Photochemical models by Moses et al. (201 1) and Line et al. 
(2010, 201 1) predict that water will be recycled in hot-Jupiter 
atmospheres, keeping H 2 O abundances close to thermochem¬ 
ical equilibrium values. A low water abundance seems to oc¬ 
cur only in atmospheres with a C/0>1. 

CO 2 , although present in hot-Jupiter atmospheres and spec¬ 
troscopically important, is not a major constituent, and it be¬ 
comes even less abundant when C/0>1. Although photo¬ 
chemistry can greatly enhance the HCN, C 2 H 2 > C 2 H 4 
abundances (Moses et al. 2013), we also see that with C/0>1, 
they are the most abundant constituents. 

In Figure 6 , the species strongly influenced by metallicity 
are again plotted as solid lines. In general, we see, as expected 
(e.g.. Line et al. 2011, Fodders & Fegley 2002, Venot et al. 
2014), that the shapes of the vertical distributions are mostly 
preserved for all metallicities. However, the thermochemical 
mixing ratio of CO 2 , CO, H 2 O, N 2 , HS, and H 2 S vary by sev¬ 
eral orders of magnitude over the range of metallicities, while 
CH 4 and hydrocarbons change very little. 

When the metallicity changes from 1 to 50, the abundance 
of CO 2 experiences the most dramatic change. It increases 
by a factor of 1000 , confirming it as the best probe of plane¬ 
tary metallicity (Fodders & Fegley 2002, Zahnle et al. 2009). 
CO 2 abundance is the quadratic function of metallicity (Venot 
et al. 2014), while CO, H 2 O, HS, H 2 S, and N 2 abundances, 
for species that either contain one metal atom or are the ma¬ 
jor reservoirs of carbon and nitrogen, increase linearly with 
metallicity (Visscher et al. 2006). For this metallicity range, 
the CO, H 2 O, HS, H 2 S, and N, abundances change by a factor 
of 100, while NH 3 , CH 4 , C 2 H 2 , C 2 H 4 , and HCN change by a 
factor of 10 or less. 

7. COMPARISON TO OTHER METHODS 

To test the validity of our code, we performed 4 different 
tests. We compared the output of TEA with the example from 
White et al. (1958) using their thermodynamic data. We also 
compared the TEA output with the output of our TEBS (Ther¬ 
mochemical Equilibrium by Burrows & Sharp) code that im¬ 
plements the Burrows & Sharp (1999) analytical method for 
calculating the abundances of five major molecular species 
present in hot-Jupiter atmospheres (CO, CH4, H20, N2, 
NH3). As another comparison, we used the free thermochem¬ 
ical equilibrium code CEA (Chemical Equilibrium with Ap¬ 
plications, available from NASA Glenn Research Center at 
http : / /WWW . grc . nasa . gov/WWW/CEAWeb/). This 
code uses the Newton-Raphson descent method within the 
Lagrange optimization scheme to solve for chemical abun¬ 
dances. Their approach is described by Gordon & McBride 
(1994), McBride & Gordon (1996), and Zeleznik & Gordon 
(1960, 1968). The thermodynamic data included in the CEA 
code are partially from the JANAF tables (Chase 1986) that 
we used in our TEA code, but also from numerous other 
sources (e.g., Cox et al. 1982, Gurvich et al. 1989, McBride 
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Fig. 7.— Left: Comparison TEA, CEA and TEBS. TEBS is an analytic method, while CEA and TEA are numerical methods. We show the major 
spectroscopically-active species in the infrai'ed that can be produced by all three methods. We run the codes for the same range of temperatures and the pressure 
of P = 1 bar. Each species in each method is plotted with a different line style, but with the same color. The TEBS final abundances are plotted as dots, CEA as 
dashed lines, while TEA is plotted as solid lines. Right: Comparison of the TEA results with the results from CEA. CEA and TEA are both numerical methods 
that use Gibbs free energy minimization method with similar optimization scheme. We show the most plausible and most abundant spectroscopically-active 
species in the infrared expected to be present in hot-Jupiter atmospheres, that all codes can cover. In the inset plot, we show a detail (zoom-in part), pointing 
out species lines that do not overlap. The T -P profile used for this run is given in the right panel of Figure 4. Tables 2 and 3 list differences between the final 
abundances for random three P.P points chosen from each run. 


et al. 1993). Lastly, we derived CEA free energies and used 
them as input to TEA, to compare the CEA and TEA outputs. 

Our first comparison was done using the example from 
White et al. (1958). We determined the composition of the 
gaseous species arising from the combustion of a mixture of 
hydrazine, N 2 H 4 , and oxygen, O 2 , at T = 3500 K and the 
pressure of 750 psi = 51.034 atm. We used the free-energy 
functions and bj values (total number of moles of element j 
originally present in the mixture) from their Table 1. We re¬ 
produced their abundances. Table 1, with slightly higher pre¬ 
cision probably due to our use of double precision. 


TABLE 1 

Comparison White ET al. vs. TEA 


Species 

RT 

White et al. 

abundances 

TEA 

abundances 

Ditference 

H 

- 10.021 

0.040668 

0.04065477 

-0.00001323 

H2 

-21.096 

0.147730 

0.14771009 

-0.00001991 

H 2 O 

-37.986 

0.783153 

0.78318741 

0.00003441 

N 

-9.846 

0.001414 

0.00141385 

-0.00000015 

N 2 

-28.653 

0.485247 

0.48524791 

0.00000091 

NH 

-18.918 

0.000693 

0.00069312 

0.00000012 

NO 

-28.032 

0.027399 

0.02739720 

-0.00000180 

0 

-14.640 

0.017947 

0.01794123 

-0.00000577 

O 2 

-30.594 

0.037314 

0.03730853 

-0.00000547 

OH 

-26.111 

0.096872 

0.09685710 

0.00001490 


Eigure 7, left panel, shows the CEA, TEA, and TEBS runs 
for the temperatures between 600 and 3000 K, pressure of 1 
bar, and solar abundances. The runs were performed with the 
input and output species that all codes contain (H, C, O, N, H 2 , 
CO, CH 4 , H 2 O, N 2 , NH 3 ). We also run the comparison just be¬ 
tween CEA and TEA, Eigure 7, right panel, for the WASP-43b 
model atmosphere that we described in Section 6 . We used 
the pressure and temperature profile shown in Eigure 4, right 
panel, and solar elemental abundances. The temperatures and 
pressures range from 958.48 to 1811.89 K and 1.5x10'^ to 


3.1623x10* bar, respectively. We included the same species 
as in Section 6 with the exclusion of the C 2 H 2 and HS species, 
because CEA does not carry the thermodynamical parameters 
for them. 

In the left panel of Eigure 7, we see that for the most species 
and temperatures CEA and TEA lines overlap (CEA result 
is plotted in dashed and TEA in solid lines). However, CH 4 
species abundances above T ~1700 K do not overlap. TEBS 
colored dots do not overplot either CEA or TEA curves, but 
follow them closely. This method is derived for only five ma¬ 
jor molecular species and is based on a few simple analytic 
expressions. 

In Eigure 7, right panel, we again see that most species over¬ 
lap, except HCN and CH 4 . The HCN curves (for CEA and 
TEA runs) differ for the full temperature range (see the inset 
figure), while, as before, CH 4 curve differs slightly only for 
pressures above ~0.1 bar and temperatures above ^1700 K 
(see Eigure 4 for the T-P profile used for this run). 

The differences seen in Eigure 7 come from the different 
sources of thermodynamic data used for CEA and TEA (see 
Tables 2 and 3). When the CEA thermodynamic data are used 
as input to TEA, all species final abundances match, see Eig¬ 
ure 8 . Section 7.1, below, elaborates on this and investigate 
the difference in free energy input values used for CEA and 
TEA. 

7.1. Comparison of free energy values in CEA and TEA 

The thermodynamic data used for CEA are in the form of 
polynomial coefficients, and are listed in the termo . inp file 
provided with the CEA code. The format of this library is ex¬ 
plained in Appendix A of McBride & Gordon (1996). Eor 
each species, the file lists, among other data, the reference 
sources of the thermodynamic data, the values of the standard 
enthalpy of formation, at the reference temperature 

of 298.15 K and pressure of 1 bar, and coefficients of specific 
heat, C^, with integration constants for enthalpy, //", and en¬ 
tropy, S”, for temperature intervals of 200 to 1000 K, 1000 to 
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Fig. 8. — Left: Comparison TEA and CEA using CEA thermodynamic data provided in their thermo. inp file. The comparison is done for the same 
conditions as in Figure 7. Tables 2 and 3 list differences between the final abundances for random three T, P points chosen from each run. 


TABLE 2 

Differences CEA v.s. TEA, Figures 7 and 8, Left Panels 


Pressure 

(bai-) 

Temp 

(K) 

CO 

CH4 

Species 

H20 

N2 

NH3 

CEA free energies 

l.OOOOe+00 

2500.00 

-33.80559930 

-34.91655970 

-40.12912409 

-27.71757996 

-32.70878374 

l.OOOOe+00 

2700.00 

-33.69182758 

-35.33414843 

-39.64186320 

-27.99507610 

-33.05270703 

l.OOOOe+00 

2900.00 

-33.61649725 

-35.76252669 

-39.25604592 

-28.25692510 

-33.39666924 

TEA (JANAF) free energies 

l.OOOOe+00 

2500.00 

-33.80793700 

-34.70780992 

-40.12426098 

-27.72037451 

-32.73695542 

l.OOOOe+00 

2700.00 

-33.69214791 

-35.08662806 

-39.63255004 

-27.99496246 

-33.08302621 

l.OOOOe+00 

2900.00 

-33.61466712 

-35.47533692 

-39.24191944 

-28.25439783 

-33.42896561 

CEA final abundances 

l.OOOOe+00 

2500.00 

5.3129e-04 

4.2666e-09 

4.3546e-04 

6.6686e-05 

8.2252e-08 

l.OOOOe+00 

2700.00 

5.2311e-04 

1.8387e-09 

4.2855e-04 

6.5661e-05 

6.4332e-08 

l.OOOOe+00 

2900.00 

5.0876e-04 

8.2340e-10 

4.1586e-04 

6.3844e-05 

4.9466e-08 

TEA final abundances using CEA free energies 

l.OOOOe+00 

2500.00 

5.3129e-04 

4.2665e-09 

4.3547e-04 

6.6685e-05 

8.225le-08 

l.OOOOe+00 

2700.00 

5.2311e-04 

1.8387e-09 

4.2856e-04 

6.5661e-05 

6.4332e-08 

l.OOOOe+00 

2900.00 

5.0876e-04 

8.2339e-10 

4.1586e-04 

6.3844e-05 

4.9466e-08 

TEA final abundances using JANAF free energies 

l.OOOOe+00 

2500.00 

5.3129e-04 

3.3976e-09 

4.3547e-04 

6.6685e-05 

8.3987e-08 

l.OOOOe+00 

2700.00 

5.2312e-04 

L4194e-09 

4.2856e-04 

6.5661e-05 

6.6260e-08 

l.OOOOe+00 

2900.00 

5.0878e-04 

6.1471e-10 

4.1586e-04 

6.3845e-05 

5.1339e-08 


6000 K, and 6000 to 20000 K. 

The JANAF tables list the reference sources of their ther¬ 
modynamic parameters in Chase (1986). The data are also 
available at http : //kinetics . nist. gov/ janaf/. 

The difference in thermodynamic parameters between CEA 
and TEA is noticeable even in their A/Z/jgg values. The 
source of standard enthalpies of formation in CEA for, e.g., 
HCN and CH 4 is Gurvich (1991), page 226 and 36, re¬ 
spectively, and their respective values are 133.08 and -74.60 
kJ/mol. The source of standard enthalpies of formation in 
the JANAF tables is listed in Chase (1986) on page 600 and 
615, respectively, and their respective values are 135.14 and 
-74.873 kJ/mol. 

TEA uses JANAF tables to calculate the values of free ener¬ 
gies for each species following Equation 10. To calculate the 
values of free energies used in CEA, we started from Chap¬ 
ter 4 in Gordon & McBride (1994). Our goal is to plug CEA 


free energies into TEA and test whether TEA will produce the 
same final abundances as CEA does. 

As explained in Section 4.2, the thermodynamic functions 
specific heat, enthalpy, and entropy as function of tempera¬ 
tures are given as; 

C" 

(35) 

HO _Jc;dT 
RT RT ’ 

(36) 

R J RT 

(37) 


These functions are given in a form of seven polynomial coef- 
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TABLE 3 

Differences CEA vs. TEA, Eigures 7 and 8, Right Panels 


Pressure 

(bar) 

Temp 

(K) 

CO 

C02 

CH4 

Species 

H20 

HCN 

NH3 

H2S 

CEA free energies 

3.8019e-01 

1719.64 

-34.9307244 

-58.4124145 

-33.595356 

-43.7404858 

-19.8124391 

-31.4721776 

-30.6054338 

1.6596e+00 

1805.28 

-34.7237371 

-57.3627896 

-33.695587 

-43.1403205 

-20.4938031 

-31.5886558 

-30.7577826 

2.1878e+01 

1810.15 

-34.7128505 

-57.3065771 

-33.701803 

-43.1083097 

-20.5310828 

-31.5955222 

-30.7664570 

TEA (JANAF) free energies 

3.8019e-01 

1719.64 

-34.9288052 

-58.4130468 

-33.5180429 

-43.7386156 

-19.6678405 

-31.4830322 

-30.5760402 

1.6596e+00 

1805.28 

-34.7231685 

-57.3648328 

-33.6061366 

-43.1392058 

-20.3573682 

-31.6020990 

-30.7285988 

2.1878e+01 

1810.15 

-34.7123566 

-57.3086978 

-33.6116396 

-43.1072342 

-20.3950903 

-31.6091097 

-30.7372812 

CEA final abundances 

3.8019e-01 

1719.64 

4.5960e-04 

5.8035e-08 

4.922 le-08 

3.768 le-04 

5.6243e-09 

7.9716e-08 

2.2498e-05 

1.6596e+00 

1805.28 

4.5918e-04 

5.2864e-08 

4.4665e-07 

3.7724e-04 

2.4131e-08 

2.8864e-07 

2.2504e-05 

2.1878e+01 

1810.15 

4.0264e-04 

5.3052e-08 

5.6873e-05 

4.3396e-04 

2.3851e-07 

3.7082e-06 

2.2520e-05 

TEA final abundances using CEA free energies 

3.8019e-01 

1719.64 

4.5959e-04 

5.8035e-08 

4.9219e-08 

3.7682e-04 

5.6240e-09 

7.9714e-08 

2.2498e-05 

1.6596e+00 

1805.28 

4.5918e-04 

5.2865e-08 

4.4667e-07 

3.7724e-04 

2.4131e-08 

2.8864e-07 

2.2504e-05 

2.1878e+01 

1810.15 

4.0263e-04 

5.3053e-08 

5.6875e-05 

4.3396e-04 

2.3851e-07 

3.7083e-06 

2.2519e-05 

TEA final abundances using JANAF free energies 

3.8019e-01 

1719.64 

4.5959e-04 

5.8326e-08 

4.5480e-08 

3.768 le-04 

4.8604e-09 

8.0472e-08 

2.2497e-05 

1.6596e+00 

1805.28 

4.5922e-04 

5.3200e-08 

4.0512e-07 

3.7719e-04 

2.0937e-08 

2.9102e-07 

2.2504e-05 

2.1878e+01 

1810.15 

4.0694e-04 

5.3429e-08 

5.2592e-05 

4.2965e-04 

2.1124e-07 

3.7386e-06 

2.2519e-05 


ficients for specific heat, Cp/R, and two integrations constants 
{ag and ag) for enthalpy, H°/RT, and entropy, S°/R\ 

C° 

—— = CL\T ^CI2T ^+^23 + cl^T + 225 ?'^ + ci(^T^ + ( 38 ) 

R 

aiT\ 


g\T) 

RT 


1/R 


St+ — 


^29S ' I ^/^298 

T \ RT 


(44) 


To see Equations 39 and 40 inside Equation 44, we multiply 
and divide the first and second term on the right with R and 
get: 


T 

-—= — a\T ^ + a2T */ 227 ’ + fl3 + fl4—+ 225—+ ( 39 ) 

K1 L j 

T^ T* as 

Q-f. -h Qi -1-, 

4 5 r ’ 


g\T) _ S"r ^ H° //»98 , . 45 . 

RT R RT RT RT ' ^ 

In the CEA analysis paper. Section 4.1, Gordon & McBride 
(1994) state that they have arbitrary assumed //°(298.15) = 
Is.fH”{29%.\5). Adopting this assumption leads to: 


5" T~^ T^ 

— = — a\ -222?' agliiTa^Ta^ -h 

R 2 2 


(40) 


^35 —h Q.'j —h Gg . 

To derive free energies in the form that TEA uses them, we 
rewrite Equation 10 for one species as: 


g^(T) 

RT 


= l/R 


-gr-^2°98l , 


RT 


(41) 


The first term on the right side can be expressed in the follow¬ 
ing format (Chase et al. 1974, Page 3): 


A)-HI 


298 




{H°-H%) 


(42) 


Thus, we rewrite Equation 41 as: 


g\T) 

RT 


= \/R 


St + 




^/^298 

RT 


(43) 


g°(T) _ S‘^ ^ ^ 

RT R RT RT RT 


(46) 


The last two terms cancel leading to a simple expression for 
free energies: 


g°{T) _S^T 
RT R RT’ 


(47) 


The first term on the right side is Equation 40, while the sec¬ 
ond term is Equation 39; expressions with polynomial coeffi¬ 
cients that are given in the CEA thermo . inp file. 

Eollowing the last conclusion, we calculated the free en¬ 
ergies for each species of interest and used them as input to 
TEA. Figure 8 shows the comparison between CEA and TEA 
using CEA free energies. We see that all species overlap. Ta¬ 
bles 2 and 3 give the exact values of free energies used and 
the final abundances for several (T,P) points that showed the 
largest differences between CEA and TEA runs in Figure 7. It 
also lists the free energies calculated using JANAF tables and 
the final abundances produced by TEA using JANAF thermo¬ 
dynamic data. 
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As seen in Figure 8, although CEA uses Newton-Raphson 
and TEA the Lagrangian method of steepest descent, both ap¬ 
proaches, using the same inputs (free energies), find the same 
final abundances. Table 2, (groups CEA final abundances and 
TEA final abundances using CEA free energies), shows val¬ 
ues identical for most species between the two tests. A few 
cases show that abundance ratios are inconsistent at the 10'^ 
level. Table 3 displays the same trend. The differences in the 
fifth decimal place may indicate that, somewhere in CEA, a 
calculation is carried out in 32-bit precision, possibly due to 
a literal single-precision number in the source code. Python 
floating literals are in 64-bit precision by default. 

8 . REPRODUCIBLE RESEARCH LICENSE 

Reproducing a lengthy computation, such as that imple¬ 
mented in TEA, can be prohibitively time consuming (Stod- 
den 2009). We have released TEA under an open-source li¬ 
cense, but this is not enough, as even the most stringent of 
those licenses (e.g., the GNU General Public License) does 
not require disclosure of modifications if the researcher does 
not distribute the code. So that the process of science can pro¬ 
ceed efficiently, there are several terms in our license to en¬ 
sure reproducibility of all TEA results, including those from 
derivative codes. A key term requires that any reviewed sci¬ 
entific publication using TEA or a derived code must pub¬ 
lish that code, the code output used in the paper (such as 
data in tables and figures, and data summarized in the text), 
and all the information used to initialize the code to produce 
those outputs in a reproducible research compendium (RRC). 
The RRC must be published with the paper, preferably as an 
electronic supplement, or else in a permanent, free-of-charge, 
public internet archive, such as github.com. A permanent link 
to the archive must be published in the paper, and the archive 
must never be closed, altered, or charged for. Details and 
examples of how to do this appear in the license and docu¬ 
ments accompanying the code, along with additional discus¬ 
sion. The RRC for this paper, including the TEA package and 
documentation, is included as an electronic supplement, and 
is also available via https : //github . com/dzesmin/ 
RRC-BlecicEtal-2015a-ApJS-TEA/. 

9. CONCLUSIONS 

We have developed an open-source Thermochemical Equi¬ 
librium Abundances code for gaseous molecular species. 
Given elemental abundances and one or more temperature- 
pressure pairs, TEA produces final mixing fractions using the 
Gibbs-free-energy minimization method with an iterative La- 


Agundez, M., Parmentier, V., Venot, O., Hersant, F., & Selsis, F. 2014, 

A&A, 564, A73, ADS, 1403.0121 

Agundez, M., Venot, 0., Iro, N., Selsis, F., Flersant, F, Flebrard, E., & 
Dobrijevic, M. 2012, A&A, 548, A73, ADS, 1210.6627 
Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 
481, ADS, 0909.0948 
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Fligh Temperature Systems: Proceedings of the First Conference, Los 
Angeles California 2-5 November 1959 (Butterworths) 

Belford, R. L., & Strehlow, R. A. 1969, Annual Review of Physical 
Chemistry, 20, 247 
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grangian optimization scheme. 

We applied the TEA calculations to several hot-Jupiter 
T-P models, with expected results. The code is tested against 
the original method developed by White et al. (1958), the an¬ 
alytic method developed by Burrows & Sharp (1999), and the 
Newton-Raphson method implemented in the free Chemical 
Equilibrium with Applications code. Using the free energies 
listed in White et al. (1958), their example, and derived free 
energies based on the thermodynamic data provided in CEA’s 
thermo . inp file, TEA produces the same final abundances, 
but with higher precision. 

Currently, TEA is specialized for gaseous species, with the 
implementation of condensates left for future work. In opac¬ 
ity calculations at low temperatures (below 1000 K), the in¬ 
clusion of condensates is necessary as it reduces the gas phase 
contribution to opacity (e.g.. Sharp & Huebner 1990, Lodders 
& Fegley 2002, Burrows & Sharp 1999). 

The thermochemical equilibrium abundances obtained 
with TEA can be used in all static atmospheres, at¬ 
mospheres with vertical transport and temperatures above 
1200 K (except when ions are present), and as a start¬ 
ing point in models of gaseous chemical kinetics and 
abundance retrievals run on spectroscopic data. TEA 
is currently used to initialize the atmospheric retrial cal¬ 
culations in the open-source BART project (available at 
https : //github . com/ joeharr4/BART). 

TEA is written in a modular way using the Python pro¬ 
gramming language. It is documented (the Start Guide, 
the User Manual, the Code Document, and this theory pa¬ 
per are provided with the code), actively maintained, and 
available to the community via the open-source develop¬ 
ment sites https://github.com/dzesmin/TEA and 
https://github.com/dzesmin/TEA-Examples. 
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